clear all;
 
%make 1D histogram for scatter density plot
y_range_b = 0; % the bottom limit for intensity
y_range_t = 8000; % the top limit for intensity
hty = [0:250:y_range_t]; %intensity ticks for 1D and 2D histogram
htx = [0 0.4];%count ticks for 1D histogram
    
[filename pathname] = uigetfile('.csv', 'select result from thunderstorm processing','multiselect','on');
for cc=1:length(filename)
    fid = fopen([pathname filename{cc}]);
    delimiter = ',';
    startRow = 2;
    formatSpec = '%f%f%f%f%f%f%f%f%f%[^\n\r]';
    ascii_data = textscan(fid, formatSpec, 'Delimiter', delimiter, 'EmptyValue' ,NaN,'HeaderLines' ,startRow-1, 'ReturnOnError', false);
    fclose(fid);
    data = [];
    Frame = [];
    IntensityS = [];
    table = [];
    x1 = [];
    y1 = [];
    data = [ascii_data{1:end-1}]; % read the .csv file 
    
    Frame = data(:,2); % Frame
    IntensityS = data(:,6); % Intensity
    [pathstr,name,ext] = fileparts(filename{cc});
    
    
    figure;
    ax1 = axes('Position',[0.3 0.3 0.3 0.2]);
    ax.ActivePositionProperty = 'outerposition';
    h = histogram(IntensityS,hty,'normalization','probability');
    TotalValue = sum(h.Values);
    ylim(htx);
    xlabel('Intensity (photons)','FontSize',12);
    ylabel('Prob','FontSize',12);
    hold on
    
    % fit a curve
    y1 = transpose(h.Values);
for ii=1:length(h.BinEdges)-1
x1(1,ii) = ((h.BinEdges(ii+1)-h.BinEdges(ii))/2)+h.BinEdges(ii);
end
x1 = x1.';

%options = fitoptions('gauss4','Upper', [Inf 700 Inf Inf 3000 Inf Inf 6000 Inf Inf 8000 Inf],'Lower', [0 500 0 0 1300 0 0 4000 0 0 6000 0]);
options = fitoptions('gauss4','Lower', [0 500 0 0 1300 0 0 2000 0 0 4000 0],'Upper',[Inf 700 Inf Inf 1500 Inf Inf 3000 1000 Inf 6000 1500]);
[gfit, gof] = fit(x1,y1,'gauss4',options);

% plot a total curve
x1 = (0:1:y_range_t);
plot(x1,gfit(x1),'LineWidth',2);
hold off
print(gcf,[[name],'fitting_curve-1.png'],'-dpng','-r300'); 

% fit first curve
    figure;
    ax1 = axes('Position',[0.3 0.3 0.3 0.2]);
    ax.ActivePositionProperty = 'outerposition';
    h = histogram(IntensityS,hty,'normalization','probability');
    ylim(htx);
    xlabel('Intensity (photons)','FontSize',12);
    ylabel('Prob','FontSize',12);
    hold on

x1 = (0:1:y_range_t);
if isempty(gfit.a1) == 0
    a1 = gfit.a1;
    b1 = gfit.b1;
    c1 = gfit.c1;
    gfit1 =@(x) (a1*exp(-((x-b1)/c1).^2));
    plot(x1,gfit1(x1),'m','LineWidth',2);
end 

% fit sencond curve

if isempty(gfit.a2) == 0
    a2 = gfit.a2;
    b2 = gfit.b2;
    c2 = gfit.c2;
    gfit2 =@(x) (a2*exp(-((x-b2)/c2).^2));
    plot(x1,gfit2(x1),'b','LineWidth',2);
end 


% fit third curve
if isempty(gfit.a3) == 0
    a3 = gfit.a3;
    b3 = gfit.b3;
    c3 = gfit.c3;
    gfit3 =@(x) (a3*exp(-((x-b3)/c3).^2));
    plot(x1,gfit3(x1),'g','LineWidth',2);
end 


% fit forth curve
if isempty(gfit.a4) == 0
    a4 = gfit.a4;
    b4 = gfit.b4;
    c4 = gfit.c4;
    gfit4 =@(x) (a4*exp(-((x-b4)/c4).^2));
    plot(x1,gfit4(x1),'c','LineWidth',2);
end 
xlim([0 6000]);
hold off
print(gcf,[[name],'fitting_curve-2.png'],'-dpng','-r300'); 


% Calculate integral area under gaussain curve

Area1 = integral(gfit1,y_range_b,y_range_t);
Area2 = integral(gfit2,y_range_b,y_range_t);
Area3 = integral(gfit3,y_range_b,y_range_t);
Area4 = integral(gfit4,y_range_b,y_range_t);
AreaTotal = TotalValue*250;
table(1,1) = a1;
table(1,2) = b1;
table(1,3) = c1;
table(1,4) = a2;
table(1,5) = b2;
table(1,6) = c2;
table(1,7) = a3;
table(1,8) = b3;
table(1,9) = c3;
table(1,10) = a4;
table(1,11) = b4;
table(1,12) = c4;
table(1,13) = Area1;
table(1,14) = Area2;
table(1,15) = Area3;
table(1,16) = Area4;
table(1,17) = TotalValue;
table(1,18) = AreaTotal;
table(1,19) = Area1/AreaTotal;
table(1,20) = Area2/AreaTotal;
table(1,21) = Area3/AreaTotal;
table(1,22) = Area4/AreaTotal;
table(1,23) = (Area3+Area4)/AreaTotal;
table(1,24) = 1-(Area1/AreaTotal)-(Area2/AreaTotal);
table(1,25) = length(Frame);
table(1,26) = max(Frame);
table(1,27) = length(Frame)/max(Frame);


% save all the parameters
tableX(cc,1) = a1;
tableX(cc,2) = b1;
tableX(cc,3) = c1;
tableX(cc,4) = a2;
tableX(cc,5) = b2;
tableX(cc,6) = c2;
tableX(cc,7) = a3;
tableX(cc,8) = b3;
tableX(cc,9) = c3;
tableX(cc,10) = a4;
tableX(cc,11) = b4;
tableX(cc,12) = c4;
tableX(cc,13) = Area1;
tableX(cc,14) = Area2;
tableX(cc,15) = Area3;
tableX(cc,16) = Area4;
tableX(cc,17) = TotalValue;
tableX(cc,18) = AreaTotal;
tableX(cc,19) = Area1/AreaTotal;
tableX(cc,20) = Area2/AreaTotal;
tableX(cc,21) = Area3/AreaTotal;
tableX(cc,22) = Area4/AreaTotal;
tableX(cc,23) = (Area3+Area4)/AreaTotal;
tableX(cc,24) = 1-(Area1/AreaTotal)-(Area2/AreaTotal);
tableX(cc,25) = length(Frame);
tableX(cc,26) = max(Frame);
tableX(cc,27) = length(Frame)/max(Frame);

fitting_table(cc).filename = [name];
fitting_table(cc).fitting_info = table;

end

save('fitting_data.mat','fitting_table');
close all

%% boostrapping for errorbar of histogram
clear;
y_range_b = 0; % the bottom limit for intensity
y_range_t = 8000; % the top limit for intensity
hty = [0:250:y_range_t]; %intensity ticks for 1D and 2D histogram

[filename pathname] = uigetfile('.csv', 'select result from thunderstorm processing','multiselect','on');
for cc=1:length(filename)
    fid = fopen([pathname filename{cc}]);
    delimiter = ',';
    startRow = 2;
    formatSpec = '%f%f%f%f%f%f%f%f%f%[^\n\r]';
    ascii_data = textscan(fid, formatSpec, 'Delimiter', delimiter, 'EmptyValue' ,NaN,'HeaderLines' ,startRow-1, 'ReturnOnError', false);
    fclose(fid);
    data = [];
    Frame = [];
    IntensityS = [];
    err = [];
    table_boot = [];
    data = [ascii_data{1:end-1}]; % read the .csv file 
    
    Frame = data(:,2); % Frame
    IntensityS = data(:,6); % Intensity
    [pathstr,name,ext] = fileparts(filename{cc});
    rng default;
    stat= bootstrp(100,@(x)x,IntensityS);
    for hh=1:size(stat,1)
        y1 = [];
        x1 = [];
        stat_his = histogram(stat(hh,:),hty,'normalization','probability','Visible','off');
        TotalValue = sum(stat_his.Values);
        y1 = transpose(stat_his.Values);
        for ii=1:length(stat_his.BinEdges)-1
            x1(1,ii) = ((stat_his.BinEdges(ii+1)-stat_his.BinEdges(ii))/2)+stat_his.BinEdges(ii);
        end
        x1 = x1.';
        options = fitoptions('gauss4','Lower', [0 500 0 0 1300 0 0 2000 0 0 4000 0],'Upper',[Inf 700 Inf Inf 1500 Inf Inf 3000 1000 Inf 6000 1500]);
        [gfit, gof] = fit(x1,y1,'gauss4',options);
        a1 = gfit.a1;
        b1 = gfit.b1;
        c1 = gfit.c1;
        gfit1 =@(x) (a1*exp(-((x-b1)/c1).^2));
        a2 = gfit.a2;
        b2 = gfit.b2;
        c2 = gfit.c2;
        gfit2 =@(x) (a2*exp(-((x-b2)/c2).^2));
        a3 = gfit.a3;
        b3 = gfit.b3;
        c3 = gfit.c3;
        gfit3 =@(x) (a3*exp(-((x-b3)/c3).^2));
        a4 = gfit.a4;
        b4 = gfit.b4;
        c4 = gfit.c4;
        gfit4 =@(x) (a4*exp(-((x-b4)/c4).^2));
        Area1 = integral(gfit1,y_range_b,y_range_t);
        Area2 = integral(gfit2,y_range_b,y_range_t);
        Area3 = integral(gfit3,y_range_b,y_range_t);
        Area4 = integral(gfit4,y_range_b,y_range_t);
        AreaTotal = TotalValue*250;
        table_boot(hh,1) = Area1/AreaTotal;
        table_boot(hh,2) = Area2/AreaTotal;
        table_boot(hh,3) = Area3/AreaTotal;
        table_boot(hh,4) = Area4/AreaTotal;
        table_boot(hh,5) = (Area3+Area4)/AreaTotal;
        table_boot(hh,6) = 1-(Area1/AreaTotal)-(Area2/AreaTotal);
        
    end
    err = std(table_boot,0,1);
    fitting_boot(cc).filename = [name];
    fitting_boot(cc).boot_info = table_boot;
    fitting_boot(cc).error_area = err;
end